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ABSTRACT 

The theory of waves and instabilities in a differentially rotating disc containing a 
poloidal magnetic field is developed within the framework of ideal magnetohydrody- 
namics. A continuous spectrum, for which the eigenfunctions are localized on indi- 
vidual magnetic surfaces, is identified but is found not to contain any instabilities 
associated with differential rotation. The normal modes of a weakly magnetized thin 
disc are studied by extending the asymptotic methods used previously to describe 
the equilibria. Waves propagate radially in the disc according to a dispersion relation 
which is determined by solving an eigenvalue problem at each radius. The dispersion 
relation for a hydrodynamic disc is re-examined and the modes are classified according 
to their behaviour in the limit of large wavenumber. The addition of a magnetic field 
introduces new, potentially unstable, modes and also breaks up the dispersion diagram 
by causing avoided crossings. The stability boundary to the magnetorotational insta- 
bility in the parameter space of polytropic equilibria is located by solving directly for 
marginally stable equilibria. For a given vertical magnetic field in the disc, bending of 
the field lines has a stabilizing effect and it is shown that stable equilibria exist which 
are capable of launching a predominantly centrifugally driven wind. 

Key words: accretion, accretion discs - hydrodynamics - instabilities - MHD - 
waves. 



1 INTRODUCTION 

There are good reasons for believing that magnetic fields are important to the physics of accretion discs. Magnetohydrodynamic 
(MHD) mechanisms provide the most convincing explanations for the anomalous transport of angular momentum that is 
required for accretion to proceed. One possibility is that angular momentum is removed from the disc by a rotating MHD 
wind (Blandford & Payne 1982). These flows have the property of collimating into jets perpendicular to the disc (Heyvaerts & 
Norman 1989), which is attractive in view of the observed association of discs and jets. A rather more convincing mechanism 
for angular momentum transport is provided by the magnetorotational instability (Velikhov 1959; Chandrasekhar 1960; Balbus 
& Hawley 1991). A Keplerian shear flow, while hydrodynamically stable, is destabilized by the presence of a magnetic field, 
provided that the magnetic field is sufficiently weak and the disc is sufficiently ionized. The instability is linear, dynamical 
and operates for toroidal as well as poloidal magnetic fields (Foglizzo & Tagger 1995; Ogilvie & Pringle 1996; Terquem & 
Papaloizou 1996). In local simulations, the non- linear development leads to turbulence, and the associated Reynolds and 
Maxwell stresses transport angular momentum radially outwards, as is required for accretion (Brandenburg et al. 1995, 1996; 
Hawley, Gammie & Balbus 1995, 1996; Stone et al. 1996). 

In an earlier work (Ogilvie 1997; hereafter, Paper I) some idealized equilibrium models of accretion discs containing 
magnetic fields were presented, with the aim of studying the magnetorotational instability in more realistic geometry. The 
model system consists of a perfectly conducting, non-self-gravitating fluid in differential rotation about a massive central 
object. The fluid contains a purely poloidal magnetic field of dipolar symmetry, which bends as it passes through the disc, 
enforcing isorotation on magnetic surfaces. This rather general model is governed by a non-linear, elliptic partial differential 
equation in two dimensions, a version of the Grad-Shafranov equation. Two classes of special solutions of this general problem 
were described. 

First, asymptotic solutions were obtained in the limit of a thin disc, using as the small parameter e a characteristic value 
of H(r)/r, where H(r) is the height of the upper surface of the disc above the equatorial plane at radius r. It was shown 
that two families of solutions resembling accretion discs exist in this limit, with different asymptotic scalings representing 
different balances of forces in the radial and vertical directions. The weakly magnetized discs are the natural generalization 
of the standard, hydrodynamic thin discs (Pringle 1981). The sound speed and Alfven velocity are both O(e) [given that the 
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Keplerian velocity is O(l)], while the fractional deviation from Keplerian rotation, caused by the radial Lorentz force, is 0(e). 
The strongly magnetized discs are not directly related to hydrodynamic thin discs but resemble models of solar prominences 
in which a sheet of matter is supported against gravity by a bending magnetic field (Kippenhahn & Schliiter 1957). The sound 
speed and Alfven velocity are both (^(e 1 / 2 ), while the fractional deviation from Keplerian rotation is 0(1). Previous approaches 
based on a vertical integration of the equations (e.g. Heyvaerts & Priest 1989) described only the strongly magnetized discs. 
It is, however, the weakly magnetized discs that are capable of launching a predominantly centrifugally driven wind if the 
magnetic field lines at the surface of the disc are inclined to the vertical by an angle greater than 7r/6. 

Secondly, solutions were obtained by assuming self-similarity in the spherical radial coordinate, as is often done in the 
analysis of winds and jets from accretion discs. This is a convenient method of studying thick equilibria and also of verifying 
the asymptotic results for thin discs. 

Although the accretion flow, the resistivity of the fluid, and any toroidal magnetic field are neglected in this model, it is 
expected that these additional terms constitute only small perturbations to the internal equilibrium of the disc. It was shown 
in Paper I that a model of a wind-driven accretion disc can be built up by superimposing these additional features on the 
solution for a weakly magnetized thin disc without disturbing the equilibrium at leading order in e. 

The primary purpose of this paper is to determine the spectrum of waves and instabilities in a weakly magnetized thin 
disc, within the framework of ideal MHD, by extending the asymptotic methods used to describe the equilibria. In particular, 
the stability of thin discs to the magnetorotational instability is to be studied. There may, of course, be other types of 
instability in magnetized accretion discs. A non-axisymmetric interchange instability (Spruit, Stehle & Papaloizou 1995) or 
bending instability (Agapitou, Papaloizou & Terquem 1997) may be present if the magnetic field provides a significant amount 
of support against gravity in the radial direction. A poloidal magnetic field is also expected to alter the criterion for convective 
instability (Moss & Tayler 1969). Finally, a global, non-axisymmetric instability, related to the Papaloizou-Pringle instability 
(Papaloizou & Pringle 1984) may be anticipated (Curry & Pudritz 1996). However, not all of these instabilities will feature 
in the analysis that follows, either because they are not present in weakly magnetized discs, or because their growth rates are 
too small to be detected. 

The energy principle of Bernstein et al. (1958) has played an important role in the analysis of the stability of magnetostatic 
equilibria relevant to astrophysics. Moss & Tayler (1969) examined the case of an axisymmetric, non-rotating star containing 
a poloidal magnetic field, and demonstrated that, if self-gravitation is neglected, the most unstable (or least stable) modes 
are those for which the azimuthal wavenumber m tends to infinity, and that the problem reduces to the consideration of a 
system of ordinary differential equations on each magnetic field line separately. Tayler (1973) considered the stability of a 
non-rotating star with a purely toroidal magnetic field, and found that the problem reduces to a system of algebraic equations 
at each separate point in the meridional plane. In each case, the most important displacements are those localized on a single 
magnetic field line or magnetic surface, and the equations obtained are closely related to those governing the continuous 
spectrum, which could be derived directly using the methods outlined in Section 2 below. 

In an accretion disc the differential rotation is the most important feature of the dynamics and cannot be ignored or 
approximated by uniform rotation. The energy principle of Bernstein et al. (1958), which applies only to static equilibria, 
cannot be used. Nevertheless, Papaloizou & Szuszkiewicz (1992) showed that, for a differentially rotating, non-self-gravitating 
fluid containing a poloidal magnetic field, a generalization of the energy principle exists in the form of a variational principle 
for the frequency eigenvalues of axisymmetric normal modes. Using the variational principle, they deduced some sufficient 
conditions for stability to axisymmetric perturbations. In this case, however, the problem does not reduce to a consideration of 
each magnetic surface separately. This suggests that the most unstable (or least stable) part of the spectrum of axisymmetric 
modes is discrete rather than continuous. This is consistent with the fact that the axisymmetric magnetorotational instability, 
although often described as a local instability, requires a relatively long wavelength in the direction perpendicular to the 
magnetic field, rather than being localized on a single magnetic surface. Nevertheless, it is possible to analyse the continuous 
spectrum directly without reference to an energy principle. 

This paper is concerned principally with the weakly magnetized thin discs and will involve an extension of the asymptotic 
methods used in Paper I. Primary consideration is given to axisymmetric modes, but non-axisymmetric perturbations are 
also discussed briefly. The continuous spectrum requires a separate analysis which can be made for a much more general 
equilibrium, and which is not restricted to axisymmetric modes. This is presented in Section 2. In Section 3 the equations and 
boundary conditions for axisymmetric modes in a weakly magnetized thin disc are derived. Some general properties of these 
equations are discussed in Section 4 and the numerical method of solution is described. In Section 5 some analytical results 
are obtained, principally for non-magnetized discs, which allow the modes to be enumerated and classified meaningfully in 
certain circumstances. For weakly magnetized discs, the most important result of this paper is the stability boundary in the 
parameter space of polytropic equilibria, which is obtained in Section 6. Non-axisymmetric modes are discussed briefly in 
Section 7, and a concluding discussion is given in Section 8. 

Throughout this paper, physical quantities are written in SI units with the permeability of free space, fio, omitted for 
convenience. The coordinates used are cylindrical polar coordinates (r, 4>, z) and the magnetic flux coordinates (tp, cj>, \) defined 



© 1994 RAS, MNRAS 000, 000-000 



Waves and instabilities in a differentially rotating disc 3 



in Paper I. The magnetic flux coordinates form a right-handed orthogonal coordinate system such that the magnetic field is 
B = W x V</> = B(ip,x) e x , where (p is the usual azimuthal angular coordinate. The Jacobian of the coordinate system is 

j = i/|vvuv<£UVx|. 



2 THE CONTINUOUS SPECTRUM 

The linearized equation governing the Lagrangian displacement £(r,t) corresponding to a small departure from any state of 
ideal MHD may be written 

p^f = -v*n - (v-€)vn - £-vvn - P |-vv$ + s v [sv£ - (v-£)B] , (2.1) 

where fl is the total pressure and 

SU = -(tp + ±B 2 )V-£ - £ ■ Vn + B ■ (B-V£) (2.2) 

is its linearized Eulerian perturbation. Here p, p, B and 7 are the density, pressure, magnetic field and adiabatic exponent, 
respectively, and D/Dt is the Lagrangian time derivative. This equation is equivalent to that given by Frieman & Rotenberg 
(1960), but is more general in that the gravitational potential $ is included, and the equation holds for an arbitrary basic flow. 
The self-gravitation of the fluid is neglected here, as it was neglected throughout the construction of equilibria in Paper I. 
Although self-gravitation has, in general, a destabilizing influence on long-wavelength perturbations, its effect on continuum 
modes, which are localized on a single magnetic surface, is nil. The equilibrium state under consideration is of the general 
class described in Section 2 of Paper I. 

The continuum modes may be described using a technique which was introduced by Papaloizou & Pringle (1982) and 
developed by Lin, Papaloizou & Kley (1993) and Terquem & Papaloizou (1996). The modes take the form of wave packets 
which are infinitely localized in the i/>-direction so that they are effectively confined to a single magnetic surface. Localization 
in the ^-direction cannot be achieved because of the infinite restoring force that would result from bending the magnetic field 
lines. Using the magnetic flux coordinates (tp, <p, \) of Paper I, one considers a sequence of functions of the form 



t(r,t) =Re 



(2.3) 



where u) € (D is the frequency eigenvalue, m G Z the azimuthal wavenumber, and k £ IR a parameter which is allowed to tend 
to infinity. The function / is a unimodal 'wavelet' 1 of width w(k) centered on tp = tpo, and forms the envelope of the wave 
packet. In the limit k — > 00, the function / is to become infinitely localized at tp = tpo, but the width w(k) of the function 
should tend to zero more slowly than fc" 1 , perhaps w oc fc -1 / 2 . Then the number of oscillations under the envelope tends to 
infinity, and a derivative of £ with respect to tp corresponds at leading order to multiplication by ik. The generalized function 
defined in the limit k — ► 00 represents an infinitely localized wave packet. It is now possible to obtain the equations satisfied 
by the mode in this limit. 

It is convenient to adopt a normalization such that |£| = O(l) at tp = tpo. If one is to obtain a solution of equation 
(2.1) with a frequency eigenvalue that has a finite limit (and a finite value of m), then it is clear that the ordering must be 
such that V-£ = O(l) [rather than O(k)] and <5II = 0(k~ 1 ) [rather than O(k) or O(l)]. Therefore the fast magnetoacoustic 
wave is filtered out in this limit, as in the magneto-Boussinesq approximation (Spiegel & Weiss 1982). It also follows that the 
component is 0(k^ 1 ), so the displacement is confined within the magnetic surface (at leading order). One may write 

V-£ = A and 511 = kT x vj, (2.4) 

where A and vj are both O(l). 

In order to write equation (2.1) in magnetic flux coordinates, it is convenient to introduce the angle of inclination i(tp, \) 
of the magnetic field to the vertical, such that 

dv L Ov 

cosi = e r ■ ej, = rB-— and sini = e r ■ e x = -r^r-^— • (2.5) 
dtp JB ox 

Then the derivatives of the unit vectors are given by 

de</> = — ( -J-7 )e x dtp + cos i d<p — ( ^- ) e x dx, (2.6) 



de^ = — cos i d(p — sin i e x d(p, (2.7) 
de x = ( -^7 ) ey, dtp + sin i e<p d<p + ( J dx, (2.8) 



1 For example, f(x) = exp(— x 2 ). 
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with 

di_ _ 1 d 

dip ~ JB d X 



and 



di_ 

dx 



= -rB 



d{JB) 
dtp ' 



(2.9) 



The quantity vj appears at leading order only in the -^-component of equation (2.1), which therefore serves to define w 
in terms of £4,, £x an d A, but is otherwise unimportant. The remaining equations at leading order are 



2ipd)fi sin i £ x — pil, £^ = — cos i r_B 



an 9$ 

lfy +p lhp 



JB 



J d X J 



(2.10) 



pu) 2 £ x + 2ip<2iQ. sin i ^ — pfi 2 sin 2 i £ x 



^£5a (J-^L\ -„h-± ( 



JB d X JB d X \JB d X J JB d X V JB d X 



1 9-0 P di> J JB \dx) J d X \J d X 



BA 



and 



= -(7p+ B 2 )A - 



t x an , b acx 



(2.11) 



(2.12) 



JB d X J d X 

where ui — u> — mil is the intrinsic frequency. These are the exact equations defining the continuous spectrum. It should be 
understood that all quantities are to be evaluated on the magnetic surface ip — tpo', the equations are essentially ordinary 
differential equations with ip treated as a parameter. Once A has been eliminated, the remaining equations may be written 

J d_ 

prJ dx 



+ 2iioQ sin i £ x 



J d X 



(2.13) 



and 



w £ x ~~ 2io>f2 sin i £0 



d 



pJB d x 



v'i + v\ 



B^_d_ 
Tdx' 



(I) 



v'i + v\ 



8 



JB 2 d X 



vl 



vi + vf 



Bg x 



where v s = (-fp/p) 1 ^ 2 and va = (B 2 /p) 1 ^ 2 are the sound speed and the Alfven velocity, respectively, 
1 d$ 2 . . 

9 * = -7Bdx~ +rn smt 

is the effective gravitational acceleration parallel to the magnetic field, and 



N x =9x 



1 dlnp g x 



JB d X 



Zx. (2-14) 



(2.15) 



(2.16) 



is a quantity analogous to the square of the Brunt-Vaisala, frequency for displacements within the magnetic surface. 

Boundary conditions must be supplied for these equations. The poloidal magnetic field line may either form a closed loop 
or extend to infinity, and in general will have segments both inside and outside the disc. The exterior region is to be treated 
as a force- free medium of zero density and pressure, in which equation (2.10) reduces to 



8_ 

dx 



(rSBj,) 



d_ 

dx 



rj_d_ 

J d X 



0. 



(2.17) 



while equation (2.11) becomes vacuous. If the field line crosses the surface 5 of the disc at points X = Xi afl d X = X2, say, 
and extends to infinity, then the relevant solution of equation (2.17) is SB^ = 0, precisely as if the exterior were a vacuum. 
Since SB must be continuous on S, the boundary condition on £4, is 

^ (^) =0 at X = Xi and x = X2 . (2.18) 

If, instead, the field line is closed, then periodic boundary conditions should be applied, with SB being continuous on S. It 
can be shown that 5B X vanishes automatically on <S (at leading order), and no boundary condition on £ x is obtained. The 
explanation is that equation (2.14) is singular at the surface, where v s vanishes, and a regularity condition applies to £ x there. 

In the absence of rotation, equations (2.13) and (2.14) are two uncoupled, second-order differential equations in Sturm- 
Liouville form, describing the Alfven continuum and the cusp continuum, respectively, for the polarization for which the 
displacement is confined to the magnetic surface. These equations have been given by Poedts, Hermans & Goossens (1985), 
who also showed that the addition of a toroidal magnetic field results in a coupling between the equations. Similarly, in this 
case, the equations are coupled by rotation, if the magnetic field is not purely vertical. A second effect of rotation is to modify 
the effective gravitational acceleration g x . 
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It should be noted, however, that there is a significant difference in the continuous spectrum between the case of no 
rotation (or uniform rotation) and that of differential rotation. In the former case, one may consider functions of the form 
(2.3) and take the simultaneous limits k — > oo and m — > oo. A continuous range of polarizations is obtained by choosing 
different limiting values of the ratio k/m, and the corresponding eigenfunctions need not have displacements confined within 
the magnetic surface. This leads naturally to the conclusions of Moss & Tayler (1969). The equations given by Poedts et al. 
(1985) in the case of a purely poloidal magnetic field describe only that part of the continuous spectrum with the polarization 
for which k/m — > oo, since m is assumed finite. When the fluid is in differential rotation, however, the limit m — > oo cannot be 
taken, because an infinite differential Doppler shift would result. The presence of a toroidal magnetic field would also prevent 
the limit m — + oo from being taken. 

It is clear that the azimuthal wavenumber m appears in equations (2.13), (2.14) and (2.18) only in the combination 
u> = ui — mQ,. Moreover, w, fi and u> are all constant on the magnetic surface. It follows that the eigenfunctions are the 
same for all values of m, while the frequency eigenvalues of non-axisymmetric modes are related to those of axisymmetric 
modes simply by a Doppler shift. In the same way that Papaloizou & Szuszkiewicz (1992) derived a variational principle for 
axisymmetric modes, so it is possible to derive a variational principle for modes of arbitrary m, provided that attention is 
restricted to the continuous spectrum. One may proceed by defining {uj 2 (tp) : n £ Z + } and {u n (ip,X) '■ n £ ^ + } to be the 
ordered, real eigenvalues and normalized, real eigenfunctions of the 'Alfven operator' that appears in equation (2.13). This is 
a Sturm-Liouville problem, 



d ( r 2 du n \ 2 2 T n ,„ im 
+ to n pr Ju n = 0, (2.19) 



dx\J d\ 

subject to the boundary condition du n /dx = at x = Xi an d X — X2 (or periodic boundary conditions for a closed field line), 
and the normalization condition 

X2 

pr 2 Ju m u n dx = 5 m n, (2.20) 

i 

for each value of ip. Let {a„} and {&„} be the components of (sini£ x /r) and (£</,/ r) with respect to this set of eigenfunctions, 
such that 



sini— = ^a„it„(V>,x) and — = ^ b n u n (tp, x)- (2-21) 

n=0 



Then the solution of the inhomogeneous equation (2.13) is given by 



ui 2 -ojl{ipy 



b - = - 75 Z:,.^ (2-22) 



and exists provided that the appropriate coefficient a n vanishes should Co 2 be equal to any of the eigenvalues of the Sturm- 
Liouville problem. The spectrum of the Alfven operator has the following characteristics. The eigenvalues {lo 2 (4>)} for each 
magnetic surface form an ordered, denumerably infinite sequence of distinct, non-negative real numbers with limit +oo. The 
lowest eigenvalue is always uJq(iP) — 0, the corresponding eigenfunction uo(ift,x) being independent of x- This zero- frequency 
mode is a trivial displacement in which the entire magnetic surface is rotated through an infinitesimal angle about the z-axis. 2 
A standard method (Courant & Hilbert 1953) leads to the asymptotic expression 

for the eigenvalues in the limit n — > oo, where 

taw= rv /2 J dx = (224) 

J xi J xi VA 

is the Alfven time of the magnetic field line. 

The eigenfunction expansions (2.21) may be substituted into equation (2.14), the equation multiplied by pJ£x an d 
integrated with respect to x to yield 

L[£ x ;uj 2 ] = R[t x ], (2.25) 



2 In the absence of a magnetic field, the infinitesimal rotation of any ring of fluid about the axis is a trivial displacement, but these trivial 
displacements are easily eliminated from the analysis by discarding the root ui = 0. 
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where the two functionals 

rx2 



and 



,2 X2 , e ,2 rq \- , ,2 



flK x ] = 4n 2 |oo| 2 + 



£ 2 



_*_V- J— * 



w 2 + «i 



(2.26) 
(2.27) 



have been introduced. Note that the infinite sum on the left-hand side starts at n = 1, the term n = 0, which is independent 
of Co 2 , having been transferred to the right-hand side. If Cj 2 = (ci> 2 ) r + i(co 2 )i, then the imaginary part of equation (2.25) is 



•'xi „ = i 



4w 2 n 2 



0, 



(2.28) 



from which it follows that the eigenvalues Co 2 are real. The eigenfunctions £ x may also be taken to be real. Now consider 
L[£ x ; ill 2 } as a function of lu 2 for a given trial function £ x . Its derivative with respect to u> 2 is 



^L[£ x ;a, 2 ]=/ p|£ x | 2 Jd X + ^ 

"'xi „ = i 



4o> 2 n 2 
(cD 2 - ^ 2 ) 2 



>0, 



(2.29) 



and so the function is monotonic wherever the derivative is defined. The function has a pole, however, at each of the values 
il> 2 — oj 2 (other than n = 0) for which the coefficient a n does not vanish. The function therefore increases once through all 
real values in each of the intervals (— oo,pi), (pi,P2), (P2,P3), ■ ■ ■ , where pi, p2, P3, ■ ■ ■ are the abscissae of the poles. Equation 
(2.25) may therefore be understood as defining a multiple- valued functional & 2 [£ x ]. Moreover, the stationary values of this 
functional may be demonstrated by standard techniques to be identical to the true eigenvalues Co 2 . Corresponding to each 
eigenfunction £ x is one eigenvalue in each of the intervals (— oo,pi), (pi,p2), (P2,P3), • • • - 3 Now L[£ x ;0] = 0, which is the 
reason for transferring the term n = to JZ[£ X ]. A necessary condition for the existence of a negative eigenvalue Co 2 is that the 
functional R[£ x ] admit negative values for trial functions satisfying the boundary conditions. The variational principle ensures 
that this condition is also sufficient. Therefore a necessary and sufficient condition for instability in the continuous spectrum is 
that the functional R[£ x ] admit negative values for trial functions satisfying the boundary conditions. In particular, stability 
is assured if 



v 2 +v 2 A/ ' 



d 



JB 2 d X 



vl + vl 



> 



(2.30) 



throughout xi < X < X2, for each magnetic surface. 

The only type of instability that can be present in the continuous spectrum is a magnetoconvective (Parker) instability 
modified by rotation. The instability is due to the release of gravitational potential energy by displacements parallel to the 
magnetic field. The first term in i?[£ x ], which may be expressed in the form 



4ft 2 |a | 2 




4ft 2 sin 2 ip|£ x | 2 Jdx, 



(2.31) 



represents the stabilizing effect of rotation, but only on modes for which ao / 0. In particular, if the disc is symmetrical about 
the equatorial plane, then this stabilizing effect is zero for all modes with odd symmetry. 

Significantly, the confinement of the displacements within the magnetic surface implies that no instabilities associated 
with differential rotation are found. In particular, the magnetorotational instability does not appear in the continuous spec- 
trum. This is in contrast with the case of a purely toroidal magnetic field (Terquem & Papaloizou 1996). Neither is a 
magnetoconvective instability associated with the -^-component of gravity found. 



3 AXISYMMETRIC WAVES AND INSTABILITIES IN THIN DISCS 

The principal aim of this paper is to determine the spectrum of waves and instabilities in a magnetized thin disc of the type 
described in Paper I. Although asymptotic analysis in the parameter e underlies the mathematics, calculations are made 
to leading order only, and therefore few explicit references to e are required. Of the two families of thin discs, the weakly 
magnetized discs have the more interesting spectrum, for they are potentially unstable to the magnetorotational instability. 



3 Although these modes share the same displacement £ x , they have different displacements according to the relation (2.22). 
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Indeed, a major objective is to determine the stability boundary in the parameter space of weakly magnetized, polytropic 
discs. Strongly magnetized discs are stable to the magnetorotational instability and so are not considered in detail. 

The analysis presented here may be seen as an extension of the work of Lubow & Pringle (1993; hereafter, LP) and 
Korycansky & Pringle (1995; hereafter, KP) on the equivalent problem for hydrodynamic thin discs. Previously, Ruden, 
Papaloizou & Lin (1988) had examined convectively unstable modes in a thin disc using a similar method. This analysis also 
extends the work of Gammie & Balbus (1994), who examined the magnetorotational instability of a stratified disc containing 
a uniform magnetic field. The present analysis is concerned with a polytropic disc containing a poloidal magnetic field that 
bends, in general, as it passes through the disc. 



3.1 Review of weakly magnetized discs 

The analysis of thin discs in Paper I is based on a small parameter e which may be defined as either the maximum value, or 
a characteristic value, of H(r)/r, where z = H(r) is the location of the upper surface of the disc at radius r. The internal 
structure of the disc is resolved by introducing a stretched vertical coordinate Q = e~ z whose value at the surface of the disc 
is C, s — e _1 H. There are two aspects to obtaining a solution at leading order. First, there are ordinary differential equations 
in £ that determine the vertical equilibrium of the disc at each radius. Secondly, there is an integral relation that determines 
the global magnetic structure; however, this is not directly relevant to the analysis that follows. 

It is assumed that the pressure and density satisfy a polytropic relation on each magnetic surface. Then the leading-order 
equations for a weakly magnetized disc are 

/GM\ 1/2 

n = (^) , (3.i) 
^ = - Po nlc+ 2porn f lBr0 , (3.2) 



(3.3) 
(3.4) 



dB r0 2p o rQ £li 
~W = B7o ' 

dQ,\ 3f2o B r Q 
lK = 2rB z0 
and 

po = K Q po, (3.5) 
where the angular velocity, density, pressure and magnetic field are given by 

Q(r,z) =fi„(r) + efii(r,C) + 0(e 2 ), (3.6) 

p(r,z) = e 2s [po(r,0 + £Pi(r,0+0(e 2 )] , (3.7) 

p(r,z) = e 2s+2 [ Po (r,0 +epi(r,C) +0{e 2 )} (3.8) 
and 

B(r, z) = e s+1 [B r0 (r, + O(e)] e r + e s+1 [B z0 (r) + 0(e)] e z . (3.9) 

Here s is a free parameter which should be positive if self-gravitation is to be unimportant. 
These equations are to be solved on < (" < Cs( r )i with boundary conditions 

B r0 (r,0) = 0, po(r, (*(r)) = and B r0 (r, ( s (r)) = B r0s (r), (3.10) 

where B r o s (r) denotes the value of B r o on the upper surface, and is supposed to be known from the global magnetic structure. 
The solution in — £ s < £ < is inferred from symmetry, p , po and fli being even functions of £, while B r o is odd. 



3.2 Normal modes 

The asymptotic method used to construct the equilibria can be applied to study their normal modes. For a weakly magnetized 
disc, the angular velocity, the buoyancy frequency, and the characteristic frequencies of acoustic and Alfven waves with 
wavelengths comparable to H are all O(l), which implies that this is the appropriate scaling for the frequency eigenvalue 
of a normal mode. Only if a mode is found with eigenvalue or oo in this scaling need a different scaling be considered. In 
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£(r,t)~ReU (r,C)exp 



— iwt + ie 



deciding the form of the eigenfunctions of axisymmetric normal modes at leading order, it is natural to assume that they take 
the form 

£(r, t) ~ Re [£ (r, C) exp(-iwt)] . (3.11) 

If this is substituted into equation (2.1), there results at leading order a sixth-order differential system 4 on — £ s < £ < £ s 
at each radius separately, which, when supplemented by appropriate boundary conditions, would constitute an eigenvalue 
problem for w. The eigenvalues would be discrete, but would depend on r in general, which means that a global solution, 
connecting a range of radii, could not be obtained other than in exceptional circumstances. One is led to conclude that the 
eigenfunction must depend more strongly on r. The correct form of the eigenfunction is 

1 J k(r)Ar | , (3.12) 

where k(r) is a radial wavenumber. This is, of course, a WKB function, but it should be emphasized that no additional 
approximation or limit is involved here; this form arises simply from the fact that the eigenfunction varies on a radial length 
scale that is much shorter than that of the equilibrium disc. 

Strictly speaking, equation (3.12) represents a wave propagating in the r-direction (when ui is real), whereas a normal 
mode is formed from the superposition of two such waves propagating in opposite directions. However, the local dispersion 
relation derived below does not depend on the sign of k, and so it is sufficient to consider travelling waves of this form. When 
u) is imaginary, equation (3.12) represents an exponentially growing or decaying normal mode rather than a travelling wave. 
The scaling of £ with e is arbitrary, since this is a linear problem, but is chosen to be O(l) for convenience. For simplicity of 
notation, the subscript '0' on all quantities is omitted hereafter. 

When the expression (3.12) is substituted into equation (2.1), the following differential system is obtained at leading 
order: 

p [(lo 2 + 3Q 2 )£, r - 2iwfi£ ] = ifcSII - V(V£ r - B r A), (3.13) 
p(u% + 2iwn&) = -V 2 ^, (3.14) 

P (- 2 + fi 2 f^) = ^? " ^ 2 CA - 2>(2>*. - B Z A), (3.15) 
where 

A = ifcCr + ^r, (3.16) 

5U = -( 7 p + B 2 )A + P Q 2 a z + B r VS, r + B Z V£ Z (3.17) 
and 

V = ikB r + B z -^. (3.18) 

This system is of sixth order and requires six boundary conditions (apart from the arbitrary normalization condition that 
applies to all linear eigenvalue problems). There are infinitely many discrete eigenvalues w for each value of k, and these form 
a local dispersion relation uj(k;r) for the disc which has infinitely many branches. The result of Papaloizou & Szuszkiewicz 
(1992) that u> 2 must be real for any global normal mode is reflected in the fact, demonstrated below, that the eigenvalues w 
of the system (3.13)-(3.15) are either purely real or purely imaginary when k is itself real. As in a standard WKB problem, a 
global solution here is trapped in a 'wave region' between two points, each of which may be either a turning point (at which 
k goes to zero) or a boundary point (either the inner or outer radius of the disc), and within which k is real and non-zero. 
The variation of £ with r is determined by higher-order terms which could be calculated in principle. 

The equations for an incompressible fluid (7 — ► 00) are obtained by imposing the constraint A = on £, and making (511 
a distinct variable in its own right, rather than using the indeterminate equation (3.17). 

The boundary conditions for equations (3.13)-(3.15) are complicated because there is a free surface with an inclined 
magnetic field. In the 'corona' above the disc, where C > Cs an d p = 0, the magnetic field is uniform on the spatial scale H 
and equations (3.13)-(3.15) reduce to the two equations 

B 2 ^=k 2 B% (3.19) 



4 Identical to equations (3.13)-(3.15) below, but with k = 0. 
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and 



ikB rs + B z -^j £* = 0, (3.20) 

where £4, is the component of £ in the meridional plane perpendicular to the magnetic field. The appropriate solution of 
equation (3.19) is 

&«exp(-|*|C), (3.21) 
when k is real. The appropriate solution of equation (3.20) is 

£4,<xexp[-i(B n /B,)kC], (3.22) 
since this corresponds to a vanishing Eulerian perturbation 

SB<t, = (ikBrs + B z -^j^=0. (3.23) 

In each case the solution is such that the components of the Eulerian perturbation of the magnetic field go to zero at infinity. 
The corresponding boundary conditions for the interior solution at £ = £ s are easily obtained by applying the continuity of £ 
and its first derivatives. Thus 

~ = -\ k \(B^ r - B„^ x ) (3 ' 24) 

and 

B,^£ = -i*Br.&. (3.25) 

The third boundary condition at £ = £s is a regularity condition, which is necessary because the sound speed goes to zero at 
the surface. The equation for £ x , the component of £ parallel to the magnetic field, has a singular point there and only the 
regular solution is acceptable. 

The three remaining boundary conditions are the equivalents of the first three at the lower surface £ = — f s . Alternatively, 
these may be replaced by symmetry conditions at f = 0, since the solutions are either even or odd. It is convenient to define an 
'even' mode as one which preserves the symmetry of the equilibrium about the equatorial plane. For such a mode, the Eulerian 
perturbation of any quantity must have the same symmetry as that quantity has in the equilibrium state. This implies that £ r 
and £4, are even functions of £, while £ z is odd. An 'odd' mode is one that breaks the symmetry of the equilibrium, implying 
that £ r and ^ are odd, while £ z is even. The boundary conditions at f = are therefore 

f = f=«.=0 ,3. 26 ) 
for an even mode, and 

Zt=U = ^=0 (3.27) 

for an odd mode. Note that this convention is opposite to that used by LP and KP. 

The relation between the local dispersion relation and the spectrum of global normal modes is quite straightforward. As 
in any WKB problem, a quantization condition of the form 

e' 1 / k(r)dr = nix + S (3.28) 

applies, where n and r2 are the limits of the wave region, n is an integer and 8 is a phase constant. This implies that 
n — 0(e _1 ) for all modes with the scalings considered, and therefore the fact that it is an integer is irrelevant at this level 
of approximation. If at any radius the local dispersion relation has a root ui corresponding to some real value of k, then a 
global normal mode must exist with that frequency eigenvalue at leading order. The spectrum is dense in the asymptotic limit 
under consideration. This also means that a necessary and sufficient condition for the existence of an unstable global normal 
mode is that the local dispersion relation should have an imaginary root ui for some real value of k. This is true, at least, for 
axisymmetric modes with dynamical [i.e. O(l)] growth rates. 

For completeness, the equations for axisymmetric waves in a strongly magnetized thin disc should be mentioned. In 
such a disc, the angular velocity and the buoyancy frequency are both O(l), but the characteristic frequencies of acoustic and 
Alfven waves with wavelengths comparable to H are both 0(e -1 ^ 2 ). This implies that an appropriate scaling for the frequency 
eigenvalue of a normal mode is 0(e -1 ^ 2 ). A set of equations is then obtained, equivalent to equations (3.13)-(3.15) but with 
O set to zero. The modes are either Alfven waves with £ r = £ z = or magnetoacoustic waves with — 0; neither rotation 
nor buoyancy has any effect, and there is no instability. However, there is a mode with exactly zero frequency in this scaling, 
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having an eigenfunction corresponding to a uniform displacement. This allows a class of modes to exist that have frequency 
eigenvalues lo — O(l) and vary on a long radial length scale, as in equation (3.11). These are the bending modes, which can be 
described using a two-dimensional treatment, and which can be unstable if the magnetic field provides a significant amount 
of support against gravity (Agapitou et al. 1997). 



4 PROPERTIES OF THE LOCAL DISPERSION RELATION 
4.1 Overview 

The sixth-order eigenvalue problem defined by equations (3.13)-(3.15) and the boundary conditions is considerably more 
complicated than the hydrodynamic equivalent. An attempt is made in Section 5 below to classify the hydrodynamic modes, 
but little progress can be made when a magnetic field is included. It is useful, however, to derive a variational principle for 
the frequency eigenvalues. This not only shows that lo 2 must be real when k is real, but also helps to clarify under what 
circumstances unstable modes are expected. The method here is an adaptation of the analysis of Papaloizou & Szuszkiewicz 
(1992) and closely follows the method of Section 2. 

Before deriving the variational principle it may be noted that the dispersion relation Lo(k) has reflectional symmetry in 
both lo- and fc-axes. This can be seen from the fact that equations (3.13)-(3.15) have the symmetries (f<£,w) i— > (— £4,, — lo) 
and either (£, of, fc) i— ► (£*, —to, —k), if lo is real, or (£,u>,k) >—> (£* ,lo, —k), if lo is imaginary. When presenting the numerical 
results it is therefore sufficient to restrict attention to positive values of lo (or Lo/i) and of k. 



4.2 A variational principle 

To invert equation (3.14), which is subject to the boundary condition (3.25), one may define the ordered, real eigenvalues 
{L0 2 (r) : n G 2 + } and the real eigenfunctions {u n {r, C) : n € Z + } of the Sturm-Liouville equation 



2 d 2 u n , 



d( 2 



+ LO„pU n = 0, 



subject to the boundary conditions 
and the normalization condition 







at C, = ±Cs 



pU m Un &C, — S„ 



(4.1) 
(4.2) 
(4.3) 



for each r. This is entirely analogous to the situation in Section 2. Physically, the eigenfunctions represent Alfven waves in 
a fictitious, non-rotating disc with the same density and vertical magnetic field, but without a radial magnetic field. The 
eigenvalues {lo 2 } are the corresponding squared frequency eigenvalues. 
The equation 



V 2 U„ + LOnpUn = 0, 

subject to the boundary conditions 
Vu n — at (= ±C S 



(4.4) 



(4.5) 



is related to the Sturm-Liouville equation by a unitary transformation. It has the same eigenvalues and its eigenfunctions are 
simply 



«n(r,C) = "n(r,C)exp 



u / (J 1 ) - 



= Un(r, C)exp 



3fi(r) 



ikr 



£> 2 ^ + L0 2 p^4, 



Equation (3.14) may be written in the form 
2kjf2p£r, 

and is conveniently solved using the eigenfunction expansions 

20 1 (r,C) 

and 



ir{r, C) = ^ a„(r)u n (r, <) exp 
and 

oo 

U(r,C) = ^6n(r)u n (r,C)exp 



3Q(r) 

2»i(r,Q 
3fi(r) 



ikr 



ikr 



(4.6) 



(4.7) 



(4.8) 



(4.9) 
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The solution is 

2iwfi(r) 



6n(r) 



uj' 2 - ul{r) 



a n {r), 



(4.10) 



and exists at any particular value of r provided that the appropriate coefficient a n (r) vanishes should J 2 be equal to any of 
the eigenvalues {co 2 (r)}. One can now substitute for £^ in equation (3.13) and form the integral relation [cf. equation (2.25)] 



involving the two functionals 



/_ . < LO^ — U) n 



(4.11) 



(4.12) 



and 



R[&,£z;k] = 4fi 2 |a | 2 + [\k\\B,Z r - BriA 2 } 



21 C= 



\m\ 



+ \V£ r \ 2 + \V£ z 



■yp + B 2 

-3pfi 2 ^| 2 -« 2 c||e- 



1 



■yp + B 2 



|B r O^ + B z ^ z +pfi 2 « z | 2 



where it has been assumed that fc is real. The imaginary part of equation (4.11) is 



' p (ICI 2 + \a 2 ) dC + (<A ]T M 2 = o, 

Sfa n=l 



(4.13) 



(4.14) 



where cj 2 = (u; 2 ) r + i(<^ 2 )i, and it follows that co 2 is real. Moreover, it can be shown that equation (4.11) has a variational 
property which implies that the necessary and sufficient condition for the existence of an unstable mode at any radius is that 
there exists a trial displacement £, satisfying the boundary conditions, which makes the functional R[(, r ,(, z ;k] negative for 
some real value of k. 

The functional R[£ r , £ z ; k] can be expressed a number of different ways which are useful for some purposes. One alternative 
is to write 

=4fi 2 |a | 2 + [\k\\B z £ r - Br£,| 2 ] C 4 

\sn\ 2 



-yp + B 2 B 2 
7P 



+ -L\B Z V& - B r V£ z \ 2 



+ 



yp + B 2 I B 



B r V£ r + B z VZ z -(—) 
\1P J 



-3 P n 2 \z r \ 2 



a dp , (pfi 2 C) 2 



7P 



ICI 2 d C . 



The coefficient of |£ z | 2 in the integrand is 



Q2 dp + (Af 



P n 2 C 



1 d In p d In p B r <9-B r 
7 9C 7P 



(4.15) 



(4.16) 



If both the specific entropy and the (radial) magnetic field strength are non-decreasing functions of £ for < £ < £ s , then 
this coefficient is non-negative and it follows that only the term — 3p$l 2 |£ r | 2 in the integrand of equation (4.15) can lead 
to instability. This implies that, in a disc with a convectively stable stratification, any instability is due to the differential 
rotation. 

A second useful variation is to extend the integral in i?[£r,£ z ; k] over the interval — oo < £ < oo and write 
R[&,£ z ;k] = 4fi 2 |a„| 2 



J — c 



'" +|2*r| 2 + |^.| a 



-yp + B 2 



1 



-3pn 2 |6.| 2 -tt 2 cg|6| 2 



■yp + B 2 
dC, 



\B r V£ r + B X V£, + pfi 2 «,| 2 



(4.17) 



etc., in which case the trial displacement must be defined on — oo < £ < oo and go to zero as |£| — > oo, but need not satisfy 
any boundary conditions at £ = ±f s . 
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Finally, the version for an incompressible fluid is 

R[^,^;k]=4Q 2 \a {) \ 2 + [\k\\B z £ r - B^] 2 ]^ + J" (W \ 2 + \V£ Z \ 2 - 3pfi 2 |£ r | 2 - Q 2 ^ \^ dC, (4.18) 

where £ is subject to the constraint A = 0. 

Despite the unconventional location of the frequency eigenvalue in the variational principle, the orthogonality principle 
is straightforward. If £ x and £ 2 are two eigenfunctions at the same radius, and for the same real value of k, which correspond 
to frequency eigenvalues uj\ and u>2, then 

(lo 2 -u; 2 2 ) f p£!-£ 2 dC = 0, (4.19) 
J-<U 

which implies that the eigenfunctions are orthogonal provided that the squared eigenvalues are distinct. Note that this 
orthogonality principle involves all three components of £, not just the meridional part. It is conjectured that the eigenfunctions 
at each radius, and for each real value of k, form a complete set in the space of continuous, vector-valued functions of C on 

-Cs < C < Cs- 



4.3 Numerical method 

The sixth-order eigenvalue problem (3.13)-(3.15) is solved numerically by the following method. First, the equilibrium itself 
must be computed, as in Paper I. When B rB and B z are specified, equations (3.1)-(3.5) constitute a third-order, non-linear 
eigenvalue problem for Q\ s (the value of fii at C = Cs)- This is solved using the shooting method: a value of fii s is guessed 
and the equations are integrated from C = Cs to C = 0. The amount by which the boundary condition B r o(r, 0) = fails 
to be satisfied defines a mismatch function whose derivative with respect to Qi s is obtained by simultaneously integrating 
the equations differentiated with respect to this parameter. Newton-Raphson convergence is then applied to obtain the 
eigenvalue. A similar technique is applied to equations (3.13)-(3.15). However, the coefficients in these equations are known 
only numerically, and so to avoid the use of interpolation one must integrate equations (3.1)-(3.5) simultaneously (with the 
previously determined eigenvalue). When k is specified, equations (3.13)-(3.15) contain only one parameter, ui. However, 
since only three boundary conditions apply at C = Csi °ne must also guess two boundary values (the sixth being given by a 
normalization condition). The full problem requires shooting in a three-dimensional complex space. However, Newton-Raphson 
iteration can still be applied successfully. 

The most difficult part of the numerical method is the treatment of the singular point at C = Cs- One can re- write 
equations (3.13) and (3.15) in the form 

H - ^ (4-20) 
c*C 7P 

d% = A, 
d( 2 -yp ' 



(4.21) 



where the numerators N r and N z are regular and are linear functions of £>, C^>, Cz and their first derivatives. The condition 
that N r and N z both go to zero as fast as 7p as C — * Cs defines a regularity condition. However, to start the numerical 
integration at C = Csi these quotients must be evaluated. This is done by expanding all quantities in power series (not Taylor 
series) about the singular point, a procedure that involves much tedious algebra. 

Once a mode has been computed for some value of k, it can be followed quasi-continuously as k is varied, so tracing out 
a single branch of the dispersion relation. This is then repeated for all branches within some limited region of the dispersion 
diagram. Fortunately there are certain limits in which the modes can be enumerated using analytical methods. 

All numerical calculations are carried out in units such that GM = Ko = r = 1, and e is defined such that Cs = 1 at the 
value of r under consideration. In particular, this means that Bo is expressed in units of 

( GM fW-i) K -Wr-%-srmr-i)g/<r-i) m (4.22) 
Taking into account the factors of e, this means that the true unit of magnetic field strength is 

(GM) r/2(r-i )jft: -i/2(r-i) r -3r/2(r-i) ff r/(r-i)_ (4 23) 

Similarly, frequencies and growth rates are expressed in units of the local Keplerian angular velocity, and the radial wavenumber 
in units of H -1 . 
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5 CLASSIFICATION OF MODES 
5.1 Hydrodynamic discs revisited 

The spectrum of axisymmetric waves in a poly tropic, hydrodynamic thin disc has been described by KP. When the magnetic 
field is absent, the solution for the equilibrium is 

ho 



Po = 



( V- 1 \ ho_ 

\ r J k 



i/(r-i) 

(5-1) 



po = K 



1\ ho_ 

K 



r/(r-i) 

(5.2) 



where the enthalpy is 

ft = i^(C s 2 -C 2 )- ( 5 - 3 ) 
Equations (3.13)-(3.15) then become equivalent to those solved by KP in terms of Eulerian perturbations. In their analysis, 
which was for a convectively stable disc, they identified three distinct classes of modes: p modes (acoustic modes, due to 
compressibility and driven mainly by pressure forces) and g modes (gravity modes, due to buoyancy and driven mainly by 
gravitational forces), both with uj 2 > Q. 2 , and also r modes 5 (inertial modes, due to rotation and driven mainly by inertial 
forces), with J 2 < Q 2 . Examples of the eigenf unctions for these modes can be found in KP. They found no modes analogous 
to the f mode (the fundamental, surface gravity mode) of stellar oscillations, but reported certain peculiarities concerning the 
ordering of the eigenvalues. The aim of this subsection is to clarify the situation and demonstrate that there are in fact two f 
modes. 

An example of the local dispersion relation for a hydrodynamic disc is shown in Fig. 1. The disc is convectively stable, 
with F = 4/3 and 7 = 5/3, and has a spectrum that is qualitatively similar to that calculated by KP using slightly different 
parameters. The different classes of modes may be distinguished (in a convectively stable disc) by the following characteristics: 

(i) f, p and g modes have uj 2 > Q 2 , while r modes have uj 2 < Q. 2 ; 

(ii) in the limit 7 — > T, in which the disc becomes adiabatically stratified, the g modes all collapse to uj 2 — Q, 2 , while all other 
modes have non-trivial limits; 

(iii) in the limit 7 — > 00, in which the fluid becomes incompressible, the p modes have lo 2 —* 00, while all other modes have 
non-trivial limits. 

The f modes are the only modes with uj 2 > fl 2 to survive both limits (ii) and (iii). There are two of them, one of each parity. 

The classification of the modes is particularly clear in the limit fc£ s — > 00. All modes with uj 2 > ft 2 are then trapped near 
the surfaces of the disc. As explained by KP, the loss of contact between the two surfaces implies that the frequency eigenvalues 
of adjacent even and odd modes must coalesce so that a single surface wave of neither symmetry can exist. The WKB method 
used by KP gives the asymptotic form of the dispersion relation for large vertical mode number, but is not appropriate for 
studying modes of small vertical mode number such as the f modes. Instead, one can use an asymptotic method based on the 
large parameter fc£ B . It is found numerically that all modes with u 2 > O. 2 are trapped in a layer of extent O ((fc^s) -1 ) near 
each surface of the disc, and have uj 2 — 0(k£ B ) in this limit. This is consistent with the interpretation of 

(i) f modes as surface gravity modes ('lj 2 ~ gk'), since g — 0(1) in the layer; 

(ii) p modes as acoustic modes ('to 2 ~ v 2 k 2 '), since v 2 — 

(iii) g modes as buoyancy modes ('ui 2 ~ iV 2 '), since iV 2 = 0(k£ B ) in the layer. 

The scalings imply that neither rotation nor the non-uniformity of gravity is significant to these modes in the limit k£ s — ► 00, 
and the modes therefore obey the same equations at leading order as in a static, polytropic atmosphere with uniform gravity. 
This problem has been solved by Lamb (1932) in terms of confluent hypergeometric functions. It was noted by Christensen- 
Dalsgaard (1980) that Lamb's analysis is valid in an asymptotic sense for modes of large degree t in an arbitrary stellar model. 
The same is true of modes of large k£ B in an accretion disc. To prove this, one should introduce a stretched coordinate 

Xl = fc(C s - C) (5.4) 

which is O(l) in the layer in which the modes are trapped. Equations (3.13)-(3.15) are then solved asymptotically with 

£r(r, C) = &-o(r,!Ci) + O ((Ks)- 1 ) , (5.5) 

Ur,0=U(r,x 1 ) + O((k(; s )- 1 ) (5.6) 
and 

lj 2 = \n 2 k( s -f-O(l), (5.7) 



5 Called g modes by LP. 
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Figure 1. Left: local dispersion diagram for a polytropic, hydrodynamic thin disc. Right: an expanded view of the lower-frequency 
branches. The disc is convectively stable, with L = 4/3 and 7 = 5/3. The frequency eigenvalues of various branches of modes, in units 
of the local angular velocity, are plotted against the radial wavenumber, in units of H —1 . The two f modes, and the first five p, g and r 
modes of each parity are shown. In order of increasing frequency these are r°, r^, rj, r|, rij, rg, r^, r|, rfj, r|, gg, gg, g|, g£, gj, g|, g|, 
§2' Sl> Si > f°> f°> Pl> Pi 1 Pl> P2' P31 P§j pIj P4' P|> pIj although even and odd f, p and g modes coalesce for sufficiently large k. 



while £<p = O ((fcCs) 1 ^ 2 ). The quantity A is the dimensionless eigenvalue of the leading-order equations, which, following 
Lamb (1932), can be combined to obtain 



xi- s ^+c^ + [2(n-l)+c-x 1 ]A o = 0, ( 5 - 8 ) 



9a;f dxi 
where 

A = i^o - (5.9) 
0x1 

is the leading-order part of V-£, while c and n are given by 

c=^-I (5.10) 
and 

LA 2 - 7 [2(r - l) n + 1] A + (7 - T) = 0. (5.11) 
The solution regular at xi = is (e.g. Erdelyi et al. 1953a) 

A oce" xl iFi(l-n;c;2a:i). (5.12) 

It is exponentially small as xi — > +00, thereby justifying the layer analysis, if and only if n is a positive integer, in which case 
the confluent hypergeometric function is proportional to a generalized Laguerre polynomial of degree n — 1. Equation (5.11) 
gives the corresponding eigenvalues A* as the roots of a quadratic equation. Provided that 7 > L > 1, these roots are real 
and ordered such that 

. .. < A2 < Ar < K A+ < \t < (5.13) 
with 

A+~ 27(r " 1) n and A" ~ 7 ~ F 1 (5.14) 
L n 27(L - 1) n y ' 

as n — > 00. The modes corresponding to A+ and A~ are called p n and g„ respectively in stellar oscillations. There is an 
additional mode which corresponds to the trivial solution A = of equation (5.8). This has 

£r0 oc £ z0 oc e~ xi (5.15) 
and A = Ao = 1. It is a surface gravity mode and corresponds to the f mode of stellar oscillations. 
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Somewhat surprisingly, the p modes and g modes are equally compressive, sharing the same function Ao, and only the f 
mode has V-£ = at leading order. Nevertheless, the g modes do survive the incompressible limit. When 7 — > 00, 

A+ -> 00 while A" -> 2( , p _ ^ - 1 ■ (5.16) 

Conversely, when 7 — > T, 

A+ -f 2(r - l)n + 1 while A~ ^ 0. (5.17) 

The notation 'f, 'p', 'g', etc., for stellar oscillations can be readily applied to thin accretion discs provided that an 
additional distinction is made between even and odd modes. This distinction is lost in the limit fc£ s — * 00, as described above. 
For general values of k, one should refer to the modes as f°, f°, p° , p°, g^j and g°, where 11 € IN is a vertical mode number 
and the superfix denotes the parity of the mode. 6 

The r modes behave quite differently in the limit k( B — > 00. In a convectively stable disc, the r modes become trapped in 
a layer of extent O ((feCs) _1/ ^ 2 ) centred on the equatorial plane, and have u> 2 = O ((fc£s) -1 ) . The asymptotic analysis proceeds 
in a similar way, with a stretched coordinate 

X2 = (fcC s ) V2 (C/Cs) (5.18) 
and scalings 

Zr(r,0~(K s )- 1/2 i;ro(r,X2), (5.19) 

e*(r,C)~£*o(r,X2), (5.20) 

e«(r,C)~e.o(r,X2) (5.21) 
and 

Lo 2 ~ (kC.^tfX. (5.22) 
The leading-order equations are 

i?H, + ^ = 0, (5.23) 

0x2 

&o = -2i\- 1/2 U (5.24) 
and 



d 2 U 



2(7 -r) 2 

T x 2 



7 (r - 1) 

A bounded solution exists if and only if 



6,0 = 0. (5.25) 



x = {2N + 1) \IW^y (5 - 26) 

where N is a non-negative integer; the solution is then 

£z0 cc cxp( — ± a 2 xI)Hn (ax 2), (5-27) 
where 

2(7-r)i 1/4 



„7(r-i)_ • 

and Hm is an Hermite polynomial (e.g. Erdelyi et al. 1953b). If iV is odd, the mode is even and may be called r° , where 
N = 2n — 1. If N is even, the mode is odd and may be called r°, where N = 2n — 2. The r modes have V-£ = at leading 
order, and their limiting behaviour depends as much on buoyancy as on inertial forces. In the marginally stable case 7 = T, 
the r modes have lj 2 = O ((fcCs)~ 2 ) an d the eigenfunctions are not localized. 

These asymptotic results are all verified numerically, as shown in Fig. 2. It should be emphasized that this classification 
scheme is based entirely on the behaviour of the modes for large values of k( B and may not accurately reflect the properties 
of the modes for smaller k^ s - 



6 An alternative classification scheme is possible in which and p° are renamed p2n and p2n— l, respectively, and similarly for other 
classes of modes. This is closer to the notation of KP, and has the advantage that the superfix is not required. However, the frequency 
eigenvalues do not then fall in sequence for all classes of modes. Also, to refer to fi and f2 rather than f° and f° is perhaps less helpful 
since it obscures the uniqueness of the f mode. 
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Figure 2. Asymptotic behaviour of the local dispersion relation for large radial wavenumbcr k, for the same disc as in Fig. 1. The f 
modes and p modes (top left) and the g modes (top right) have u) 2 = O(k) as k — > oo, while the r modes (bottom) have up = 0(fc _1 ). 
In each case the dotted lines indicate the asymptotic limits derived in the text. Note that the vertical axis is different in each plot. 



5.2 Magnetized discs 

The presence of a poloidal magnetic field in the disc has profound implications for the spectrum of waves and instabilities. As 
well as modifying the modes that occur in a hydrodynamic disc, the magnetic field gives rise to additional modes which are 
potentially unstable. Furthermore, the branches are no longer separated in the dispersion diagram, but undergo avoided cross- 
ings. This makes the classification of modes difficult, although there is one limit which can be investigated semi-analytically. 
This is the case k = in a disc with a purely vertical magnetic field, as was considered by Gammie & Balbus (1994). 

When the magnetic field is purely vertical, the Lorentz force vanishes and equations (5.1)-(5.3) for the equilibrium are 
valid. If also k — 0, equations (3.13)-(3.15) simplify considerably because the horizontal and vertical components of £ become 
decoupled. There exist purely horizontal modes, with 

£ r {r,0 = a r (r)u N (rX) (5.29) 
and 

^(r,C) = O0(r)ujv(r,C), 
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where un is one of the eigenfunctions of equation (4.1), and the components a r and satisfy the equation 



tu 2 + 3D. 2 - u 2 N -2iwfi " 




a r 




~o~ 


2iu>n uj 2 — ui% 




a.0 








The frequency eigenvalues are the roots of the equation 

lj 4 — (2lo 2 n + Q 2 )cu 2 + lu 2 n (lu 2 n — 3S1 2 ) = 0, 



(5.31) 



(5.32) 



2 i 1 

- <JJ N + a" 



1 ± (l + 16lu 2 n /Q 2 ) 



1/2 



(5.33) 



so that there is an unstable mode if and only if < uj 2 n < 3Q, 2 . Since the eigenvalues {lu%} are ordered, the JV = 1 mode 
is unstable if any mode is unstable. The maximum possible growth rate is the Oort parameter A = which is achieved if 

0, the eigenvalues uj 2 n — > (although not uniformly with respect to JV), and equation (5.33) 



u) 2 N 



= ifO 2 . In the limit B z 



reduces to 



fT or 0. 



(5.34) 



Accordingly, the lesser solution of equation (5.33) may be called an m mode, since it is due entirely to the magnetic field. 
Specifically, it is m^ if JV = 2n, or if JV = 2n — 1. If JV = the mode is trivial. The greater solution of equation (5.33) 
may be called r® if N — 2n, r° if N = 2n — 1, or P if N = 0. However, the r and g modes really lose their identity when a 
magnetic field is introduced. 

There are also purely vertical modes, which satisfy the equation 



d( 2 



{l + 2q)C^ + N{N + 2q)Z z = 0, 



where 



q = 



r + i 

2(r - 1) 



and 



(5.35) 



(5.36) 



(5.37) 



Appropriate solutions, regular at £ = ±f s , are obtained when N is a non-negative integer, in which case 

e»occ&(c/&), ( 5 - 38 ) 

where is a Gegenbauer polynomial (e.g. Erdelyi et al. 1953b). Equation (5.37) then leads to the dispersion relation 

r + i 



jV 



N ■ 



(5.39) 



These modes may be identified as if JV = 2n - 1, p° if JV = 2n, or f ° if JV = 0. 7 

An example of the dispersion diagram for a disc containing a purely vertical magnetic field is shown in Fig. 3. In this 
case the strength of the magnetic field is such that one unstable magnetorotational mode of each parity exists. These modes 
are unstable only for sufficiently small k, and their growth rates are greatest for k = 0. The stable modes form a complicated 
dispersion diagram because of the avoided crossings that occur between different branches. Broadly speaking, the pattern 
consists of branches which, like the f and p modes in a hydrodynamic disc, sweep upwards in the diagram, and other branches 
which are almost flat. Where these would cross over each other, a closer examination reveals that each branch is in fact 
continuous; the two branches approach and then diverge without crossing. The character of the eigenfunctions is exchanged 
smoothly in the neighbourhood of these avoided crossings. This means that a classification of the modes based on continuity 
with k — would not be meaningful. The reason for plotting even and odd modes separately is that avoided crossings occur 
only between modes of equal parity. The fact that parity is a discrete property which cannot be exchanged smoothly implies 
that modes of opposite parity cannot interact in this way, and so branches of even and odd modes cross freely over one another. 
Similar avoided crossings in the spectrum of a magnetized isothermal atmosphere (without the symmetry of reflection in a 



7 In fact the f° and r ° modes in a Keplcrian disc without a magnetic field involve both horizontal and vertical displacements even when 
k = 0. 
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Figure 3. Part of the local dispersion relation for a thin disc containing a purely vertical magnetic field. The parameters are T = 4/3, 
7 = 5/3 and B z = 0.01. Modes with even and odd symmetry are shown in the left and right panels, respectively. The frequency 
eigenvalues of various branches of modes, in units of the local angular velocity, are plotted against the radial wavenumber, in units of 
H -1 . A dashed line indicates an unstable mode for which oj is imaginary. In this case one unstable magnctorotational mode of each 
parity exists. A maximum growth rate of 0.7151C2 is achieved at k = by the mj mode. The frequency eigenvalues of the continuous 
spectrum, which are all real, are marked on the scale at the right of each panel; these are the eventual limits of the branches as k — > oo. 




02468 10 02468 10 

k k 



Figure 4. Part of the local dispersion relation for a thin disc containing a purely vertical magnetic field. The parameters are T = 4/3, 
7 = 5/3 and B z = 0.002. Modes with even and odd symmetry are shown in the left and right panels, respectively. In this case seven 
unstable magnctorotational modes of each parity exist. A maximum growth rate of 0.7495O is achieved at k = by the m| mode. The 
frequency eigenvalues of the continuous spectrum, which are all real, are marked on the scale at the right of each panel. 

horizontal plane) have been analysed by Hasan & Christensen-Dalsgaard (1992). 

When the magnetic field is weaker, more unstable modes exist but their branches fit neatly without mode interactions, 
as shown in Fig. 4. It is well known that the addition of a weak magnetic field to a hydrodynamic disc constitutes a highly 
singular perturbation if ideal MHD is assumed. This is reflected in the fact that the convergence of the frequency eigenvalues 
of m n modes to zero as B — > is not uniform with respect to n. Modes of large n clutter the dispersion diagram even when 
the magnetic field is exceedingly weak, and can be unstable with large growth rates. 

Finally, an example of the dispersion diagram for a disc containing a bending poloidal magnetic field is shown in Fig. 5. 
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Figure 5. Part of the local dispersion relation for a thin disc containing a bending poloidal magnetic field. The parameters are T = 4/3, 
7 = 5/3, B z = 0.015 and B rB = 0.005, with Q± B 0.3376. Modes with even and odd symmetry are shown in the left and right panels, 
respectively. In this case only the mj mode is unstable, achieving a maximum growth rate of 0.6019O at k = 0. The frequency eigenvalues 
of the continuous spectrum, which are all real, are marked on the scale at the right of each panel. 

The avoided crossings are wider than for a purely vertical magnetic field, which implies that the coupling between different 
branches of modes is stronger when the field lines bend. 

The behaviour of branches in the limit k — » oo is quite different from the hydrodynamic case. If any one branch is followed 
continuously, it undergoes a finite number of avoided crossings and then co approaches a finite limit. It can be shown that 
these limits are the frequency eigenvalues of the continuous spectrum discussed in Section 2. In this limit the Lagrangian 
displacement is confined entirely within the magnetic surface and the magnetorotational instability, which is associated with 
the gradient of angular velocity perpendicular to the magnetic field, is lost. 



6 STABILITY CRITERIA FOR MAGNETIZED DISCS 

The stability of a weakly magnetized thin disc to axisymmetric perturbations can be decided, in principle, by computing the 
local dispersion relation at each radius separately; the disc is unstable if and only if an imaginary frequency eigenvalue exists, 
for any real value of k, at any radius. However, a more efficient algorithm is required in practice. In the case of polytropic 
discs, the parameter space divides into stable and unstable regions separated by a curve of marginal stability. The location 
of the marginal curve depends on a further parameter which is the adiabatic exponent 7. Equilibria on the marginal curve 
possess a mode with u> = 0, but this is not sufficient to determine the curve, since all unstable equilibria also possess modes 
with u> — (cf. Figs 3-5). In principle, a marginal curve can be drawn for each m mode and for every value of k. 

Before describing the marginal curves, it is appropriate to give a more detailed account of the solutions of the equilibrium 
equations for polytropic discs. When B rB and B z are specified, equations (3.1)-(3.5) constitute a non-linear eigenvalue problem 
for f2i s . The solutions of these equations lie on a two-dimensional manifold in the three-dimensional parameter space of 
(B rB , B z , ills). Some of these solutions are physically acceptable and may be called 'regular' equilibria: the magnetic field 
lines bend only once when passing through the disc, and the density and pressure decrease monotonically from the equatorial 
plane to the surface. In the remaining solutions, which may be called 'irregular' equilibria, the field lines bend more than 
once; some of these solutions also have density inversions. In Paper I only the class of regular equilibria (described as the 
'principal branch') for the case T = 5/3 was mentioned. It is now to be shown that all irregular equilibria are unstable, so 
that the marginal curves need be drawn only on the principal branch. 

Theorem 1. All irregular equilibria possess an unstable mode of odd symmetry at k = 0. 

Proof. When k = 0, the form of the functional R[£ r , £*; k] for purely horizontal trial displacements of odd symmetry is, from 
equation (4.13), 

-3 P r2 2 |£,| 2 ) dC, (6.1) 



R[&,o-,o] 
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and the trial function must satisfy d£ r /dC, = at f = ±£ s . Given that the equilibrium is irregular, B r must have zeros, other 
than at £ = 0, at points £ = ±C n , where < £i < £2 < • • • < Ov < Cs aQ d iV > 1. Then consider the Sturm-Liouville equation 

B 2 ^f+A P2 / = (6.2) 

on — ("i < C < (1, subject to the boundary conditions 

|U0 at C = ±Ci- (6-3) 

Since the equation is symmetrical about f = 0, the eigenfunctions j/o, J/i, J/2, • • • alternate between even and odd symmetry, 
and the eigenvalues A , Ai, A2, ... form an ordered, increasing sequence of distinct, non-negative real numbers. The lowest 
eigenvalue is Ao = 0, corresponding to y a = constant. Now it follows from equations (3.1)-(3.4) for the equilibrium that fii 
satisfies the equation 

Bf^pi+3pn 2 n! = o, (6.4) 

and also the boundary conditions (6.3). Moreover, since (1 is the first zero of B r in < f < f s , B r has exactly two extrema, 
and £7i has exactly two zeros, in the domain of the Sturm-Liouville equation. It follows that f2i is the eigenfunction t/2, and 
has eigenvalue A2 = Sfl 2 . The properties of the Sturm-Liouville equation then ensure that the first odd eigenfunction y\ has 
an eigenvalue satisfying the (strict) inequalities < Ai < Zfl 2 . Now construct the trial function 

f(t\ /^(O. ICI<Ci, (R ~ 

/(f) = \ BB n(0in(Ci), ICI>Ci, (6 ' 5) 

which is continuous throughout — £ s < £ < £ s and satisfies the boundary conditions for the variational principle. After an 
integration by parts one obtains 

fl[/,0;0] = -2 / >Cl (3 P n 2 -A 1 )|/| 2 dC-2 [' ipfl 2 |/| 2 dC, (6.6) 
Jo JCi 

and, since this is clearly negative, the theorem is proved. 

In Fig. 6, the different branches of equilibria, now for the case F = 4/3, are shown in a number of suitably chosen sections 
through the parameter space. It is found that, as B z is reduced, an increasingly large number of branches of irregular equilibria 
appear as the solution manifold folds and divides into many leaves. There is a critical value of B z , approximately 0.01283, 
below which the principal branch becomes separated from the 'vertical' solution (which has B r = and flu = 0) by a stretch 
of irregular equilibria. (The vertical solution, however, is technically always regular.) As B z is reduced further, a bifurcation 
disconnects the principal branch from the vertical solution. An infinite number of further bifurcations occur as B z — ► 0. 

It is observed in the numerical calculations that, for regular equilibria, the last mode to be stabilized is the m° mode, 
and it is stabilized last at k = 0. This supports the following conjecture. 

Conjecture. If a regular equilibrium possesses an unstable mode at any value of k, then it possesses an unstable mode at 
k = 0, assuming that 7 > T. 8 

While this conjecture has not yet been proved in the most general circumstances envisaged, the following two theorems support 
it strongly. 

Theorem 2. The conjecture holds in the case of an incompressible fluid, provided there is no density inversion. 
Proof. An appropriate version of the functional -R[£ r ,£ z ; k] for an incompressible fluid is 

R[^^ z ;k]=4fl 2 \a () \ 2 + J (\V£ r \ 2 + \V£ Z \ 2 -3pfl 2 \Z r \ 2 -fL 2 (^ z \ 2 ^j d<, (6.7) 

where the trial displacement £ is defined on —00 < £ < 00 and is subject to the constraint A = 0. Given that there exists an 
unstable eigenfunction at some wavenumber k =fc 0, it follows that 

R[ir,iz;k]<0. (6.8) 

When k = 0, the constraint A = can be satisfied by using a purely horizontal trial displacement £ r = /, where the function 
f(() is defined by 

£■(0 = /(C) cxp [-(2ni/3n)iAr] , (6.9) 



8 When 7 < r there may be a (magneto-) convective instability and the result cannot be expected to hold. In fact, a condition differing 
slightly from 7 > T may be required in order to prove the conjecture. 
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Figure 6. Solution curves for weakly magnetized, polytropic thin discs with V = 4/3. The 'eigenvalue' f2i s is plotted against B rs for 
various values of B z (as indicated in each panel). All solutions in the range —0.025 < B rB < 0.025 arc shown. Regular and irregular 
equilibria are indicated by solid and dotted lines, respectively. As B z is decreased, the solution manifold folds and divides into many 
leaves, and the branch of regular equilibria becomes separated from the solution at the origin which has a purely vertical magnetic field. 



so that 



V (k) i r = ikB r i r + B z ^ = B z (^j cxp [- (2fii /3n)iAr] = (z> (0 ) /) cxp [~(2ni/3n)iAr] , 



(6.10) 
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where the suffix on the operator V indicates the value of k used. Then 



Hr\ 2 = |/| 2 and \V (k) i 



l^(0)/| 2 - 



(6.11) 



Also, the coefficients {a„} in the eigenfunction expansion (4.8) are the same when £ r = / and k = as they are when £ r = £r 
and k — k. Thus the difference 



«Kr,6;fc]--R[/,0;0] = 

is non-negative, and so 

R[f, 0;0] < R[ir,i z ;k] < o, 

which proves that an unstable mode must exist with k — 0. 

Theorem 3. The conjecture holds for a compressible fluid when the magnetic field is purely vertical. 
Proof. In this case one may take 
=4fi 2 |a | 2 



(6.12) 



(6.13) 



l<5n| ; 



yp + B 2 



+ 5 



IP 



-yp + B 2 



B l 



dC, \ 7P 



-3pfi 2 |^| 2 



where 

sn = -(tp+ b 2 )a + P n 2 ^ z + b 2 ^. 



d( 7P 



Given that there exists an unstable eigenfunction £(£) at some wavenumber k 0, it follows that 
R[£r,£,;k] <0. 



(6.14) 



(6.15) 



(6.16) 



When = 0, the purely horizontal displacement £ r = £ r should be used, and 511 then vanishes. The coefficient of |£ z | 2 in the 
integrand of equation (6.14), 

' 7 -r\ (p« 2 c) 2 



oC 7P 



-(^) 



7P 



(6.17) 



is non-negative, since it is assumed that 7 > T. Thus 
R[ir,0;0] < R[£r,£z;k] < 0, 

which proves that an unstable mode must exist with k — 0. 

If the conjecture is accepted on the basis of these lemmata and the numerical evidence, it implies that the overall stability 
boundary in the parameter space is the marginal curve for a mode (in fact, the mj mode) at k = 0, drawn on the principal 
branch. This curve has been computed directly by solving the equations for an equilibrium (with a given value of B z ) and a 
mode (with ui — and k = 0) as an eigenvalue problem, yielding a value for B rs - The stability boundaries for equilibria with 
r = 4/3 and T = 5/3 are shown in Fig. 7, using 7 = 5/3 for the adiabatic exponent. In comparing the two panels it is important 
to note that the unit of magnetic field strength depends on Y, according to equation (4.22). It is more convenient, therefore, to 
compare a dimensionless quantity such as the plasma beta [3 = 2p/B 2 . The plasma beta of the marginal equilibria, evaluated 
on the equatorial plane, is plotted against B z in Fig. 8. It is seen that all marginal equilibria have values of /3 reasonably close 
to unity. However, f3 is not generally a reliable guide to stability, especially for equilibria in which the angle of inclination of 
the magnetic field exceeds 7r/6. 

The emergence and disappearance of equilibrium solutions as the parameters are varied is not necessarily governed by 
conventional bifurcation theory. Typically, an equilibrium ceases to exist when the parameters are changed in such a way 
that would make the enthalpy become negative at some point. This would correspond to a branch point of the polytropic 
relation (3.5), and non-analytic behaviour is to be expected. In fact, the equilibrium may continue to exist in a complex- valued 
sense, even though this has no physical meaning. This means that there is little hope of using the bifurcations of the solution 
manifold to infer the stability of the solutions. Indeed, two further obstacles to such a method present themselves: first, only 
the stability to modes of even parity could be considered, since the equilibria are all symmetric; and secondly, the effect of 
the adiabatic exponent 7 on stability could not be taken into account. 

The magnetorotational modes are not affected greatly by buoyancy or compressibility. In Fig. 9 the stability boundaries 
for equilibria with Y = 4/3 and Y = 5/3 are plotted for different values of the adiabatic exponent. The unstable region reduces 
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Figure 7. Stability boundaries in the parameter spaces for weakly magnetized, polytropic thin discs with T = 4/3 (left) and T = 5/3 
(right). In each case, the dotted line indicates an angle of inclination i = n/6. The dashed line is the critical curve for the existence 
of the principal solution branch. The solid line is the marginal curve, at k = 0, for the m° mode, when 7 = 5/3. This is the overall 
stability boundary of the equilibria to axisymmetric perturbations. In the regions Ft, which extend indefinitely beyond the upper right 
of each figure, there exist equilibria which arc capable of driving a wind and are also stable to the magnetorotational instability. Note 
the different scales used in each plot. 
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Figure 8. The plasma beta of the marginal equilibria, evaluated on the equatorial plane, plotted against the vertical magnetic field 
strength in the interval in which a marginal equilibrium exists. The two panels correspond to the two panels of Fig. 7, i.e. equilibria with 
T = 4/3 (left) and with T = 5/3 (right), where 7 = 5/3 in each case. Marginal equilibria typically have values of /3 that are reasonably 
close to unity. 

slightly in size as 7 is increased from F to 00, demonstrating the mildly stabilizing effect of a sub-adiabatic stratification. The 
effect is most pronounced for equilibria in which the magnetic field bends significantly, and is zero for those with a purely 
vertical magnetic field. 



7 NON-AXISYMMETRIC MODES 

The analysis of Section 3 may be extended to include non-axisymmetric waves and instabilities of azimuthal wavenumber m. 
A distinction must be made between modes with m = 0(e _1 ) and those with m = 0(1). In the former case, the azimuthal 
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Figure 9. The effect of sub-adiabatic stratification on the magnctorotational instability. Left: stability boundaries for weakly magnetized, 
polytropic thin discs with Y = 4/3 are plotted for four different values of the adiabatic exponent 7. The four solid lines are the stability 
boundaries for 7 = 4/3, 5/3, 3 and 00 (from right to left). Right: similar results for equilibria with T = 5/3. The three solid lines are the 
stability boundaries for 7 = 5/3, 3 and 00 (from right to left). 



wavenumber is comparable to the radial and vertical wavenumbers, and the equations of Section 3 are not valid in any sense. 
These modes are expected to resemble the global non-axisymmetric instabilities of thick tori studied by Papaloizou & Pringle 
(1984). The strong differential Doppler shift means that these modes are probably localized in a region of small radial extent 
about their corotation radius, and require a boundary of the disc to be present in this small region if they are to have 
dynamical [O(l)] growth rates. 

For modes with m = O(l), however, the equations of Section 3 require very few changes. The eigenfunctions have the 
form 



£(r,t)~ReU (r,C)exp 



-iujt + \mS + 



ie 1 J k(r) dr j> 



(7.1) 



and the frequency eigenvalue u) must be replaced with the intrinsic frequency w = u) — mil. In this way LP were able to discuss 
the propagation of non-axisymmetric waves in an isothermal accretion disc. Provided that the dispersion relation yields real 
values of uj for real values of k, there is no difficulty because cu, u) and k are all real within the wave region of a mode. 

There remains the possibility that the dispersion relation yields imaginary values of to for real values of k, as a result of 
the magnetorotational instability or a convective instability. In that case one cannot simply replace to with u>, for the following 
reason. Since u> must be constant for a normal mode, it follows that the real part of u) can vanish at (at most) one radius, the 
corotation radius of the mode. Therefore w is imaginary at the corotation radius, but anywhere else it cannot be a solution 
of the dispersion relation for any real value of k. The conclusion is that, for a non-axisymmetric unstable mode, k is real only 
at the corotation radius. 

When k is not real, the variational principle of Section 4.2 is not valid. Otherwise, the only change that need be made in 
Section 3 is to replace |fc| with the function 



p(k) 



+k, Re(fc) > 0, 
-k, Re(fc) < 0, 



(7.2) 



which is an analytic function apart from a branch cut along the imaginary axis. It then follows that each branch of the 
dispersion relation is an analytic function, except on Re(fc) = 0. The derivative dui/dk is known for real values of k from the 
real dispersion relation, and is imaginary for an unstable branch. Let k — k c be the real value of the radial wavenumber at 
the corotation radius r = r c of such a mode. Then, at a neighbouring radius r, the wavenumber is (to first order) 



fc c + 



3mQ c 
2r c 



dkJc 



(r - r c ), 



(7.3) 



assuming that the derivative (duj/dk) c does not vanish. The quantity in square brackets is purely imaginary; depending on 
the sign of the various terms, Im(fc) is either an increasing function or a decreasing function of (r — r c ). If it is an increasing 



© 1994 RAS, MNRAS 000, 000-000 



Waves and instabilities in a differentially rotating disc 25 

function, then the mode will take the form of a localized Gaussian wave packet around r — r c , according to equation (7.1). If 
it is a decreasing function, the function (7.1) would increase greatly in magnitude away from r = r c and cannot satisfy the 
radial boundary conditions. 

In conclusion, non-axisymmetric unstable modes are localized around the corotation radius, and obey the real dispersion 
relation only at that point. The modes have either positive or negative values of m depending on the sign of the imaginary 
group velocity. At any given point on an unstable branch of the real dispersion relation, non-axisymmetric modes of this kind 
can be found. 

8 DISCUSSION 

This paper has addressed a number of issues relating to the spectrum of waves and instabilities in an accretion disc containing 
a poloidal magnetic field. The general analysis of the continuous spectrum in Section 2 demonstrates that the only type 
of possible instability that is truly localized on a single magnetic surface is an essentially axisymmetric magnetoconvective 
instability associated with the component of gravity parallel to the magnetic field. Unlike the case of uniform (or zero) rotation, 
the interchange instability, if present, is not manifest in the continuous spectrum. Neither is the magnetorotational instability 
or any other instability associated with differential rotation. These must instead be sought in the normal modes of the system. 

The extension of the asymptotic methods of Paper I leads to a WKB description of axisymmetric waves and instabilities 
in a weakly magnetized thin disc in terms of a local dispersion relation which generalizes earlier work by LP and KP. The 
waves propagate radially as in a slowly varying waveguide. The modes in a hydrodynamic disc can be classified as f, p, g 
and r modes according to their behaviour in the limit of large radial wavenumber. When a magnetic field is introduced, the 
dispersion diagram is complicated by a large number of avoided crossings which make a systematic classification difficult. If 
the magnetic field is relatively weak, unstable magnetorotational modes also occur for sufficiently small values of the radial 
wavenumber. The overall stability boundary in the parameter space of polytropic equilibria is obtained by computing the 
marginal curve for the first magnetorotational mode of odd symmetry. For a given value of the vertical magnetic field, the 
addition of a radial component (which makes the field lines bend) has a stabilizing influence. A sub-adiabatic stratification 
also has a mild stabilizing influence if the magnetic field is not purely vertical. 

This analysis shows for the first time that it is possible to construct stable equilibria which are capable of driving a wind. 
Indeed, increasing the magnetic field strength not only tends to stabilize the equilibria but also makes it easier to construct 
equilibria in which the magnetic field lines at the surface of the disc are inclined to the vertical at angles significantly greater 
than 7r/6. It should be emphasized, however, that the stability analysis applies only to the disc and not to any wind solution 
that may be superimposed on it. 

It is also important to note that the analysis is valid only for axisymmetric modes and non-axisymmetric modes of 
small azimuthal wavenumber [m = O(l)], and that it can only detect instabilities with dynamical [O(l)] growth rates. It is 
expected that all models are subject to global non-axisymmetric instabilities which depend strongly on the radial boundaries 
of the disc, resembling the Papaloizou-Pringle instability, but these probably have very small growth rates in a thin disc. The 
radial interchange instability (Spruit et al. 1995), if it is not entirely stabilized by the differential rotation, would also have 
a sub-dynamical growth rate in a weakly magnetized thin disc. The equilibria described here as 'stable' may therefore not 
exist as truly laminar flows but could be subject to weak turbulence or at least fluctuations. However, that would be quite 
different from equilibria that are locally unstable to the magnetorotational instability, which are bound to degenerate into 
strong MHD turbulence. 

There are many ways in which this analysis could be improved and extended. In order to study the radial propagation of 
waves, an explicit global equilibrium model of a magnetized disc should be constructed by choosing the functional forms of 
Cs(f)j V'o(f) and Ko(r), say, and solving for the equilibrium at each radius. The variation of the radial wavenumber with radius 
could then be determined for any mode by following the dispersion relation. The dependence of the amplitude of the mode 
on radius could be found either by solving the linearized equations at the next order in e, or, more simply, by appealing to a 
wave-action conservation relation (cf. LP). Other areas to be explored are the possible stabilization of a super-adiabatically 
stratified disc by a sufficiently strong magnetic field, the propagation of the m = 1 'tilt' mode in a magnetized disc, and the 
effects of self-gravitation on the equilibria and spectra of both hydrodynamic and magnetized discs. 
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